%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     Lab 3: Digital Filter Design                        %
%                   EE558L: Section 1 (F: 1200-1440)                      %
%                            Dr. Nagaraj                                  %
%                     Author: Michael Spinali                             %
%                             813488956                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Create Fresh Environment
close all
clear all
clc

%%%%%%%%%%%%%%% PART 1 %%%%%%%%%%%%%%%
% Define Constants
N = 15;                             % Length of Filter
fs = 16E3;                          % 16kHz Sampling Frequency
f = [0 2.5E3 2.7E3 8E3]./(fs/2);    % firpm freq. vector
a = [1 1 0 0];                      % firpm freq. amplitude vector

% Loop once for 15-0rder, and once for 31-0rder
while(N <32)
    
    % Create Filter
    b = firpm(N,f,a);
    
    % Plot Filter in Time Domain
    figure;
    plot(0:N,b);
    grid on;
    axis([0 N -0.1 0.35]);
    title(sprintf('h[n] (%d-order LPF)',N),'fontsize',14);
    xlabel('Time (samples)','fontsize',14);
    ylabel('Amplitude','fontsize',14);

    % Determine Freq. response of filter
    figure;
    [h,w] = freqz(b,1,1024);
    
    % Plot Magnitude
    subplot(2,1,1);
    plot(f.*(fs/2),a,'r',(w.*(fs/2))/pi,abs(h),'b')
    grid on;
    legend('Ideal','firpm Design')
    title(sprintf('Magnitude Plot (%d-order LPF)',N),'fontsize',14);
    xlabel('Frequency(Hz)','fontsize',14);
    ylabel('|H[f]|','fontsize',14);

    % Plot Phase
    subplot(2,1,2);
    plot((w.*(fs/2))/pi,angle(h),'b')
    grid on;
    title(sprintf('Phase Plot (%d-order LPF)',N),'fontsize',14);
    xlabel('Frequency(Hz)','fontsize',14);
    ylabel('<H[f]','fontsize',14);
    legend('firpm Design')
    N = N+16;
end

%%%%%%%%%%%%%%% PART 2 %%%%%%%%%%%%%%%
f1 = 1.5E3;                 % Freq 1
f2 = 6E3;                   % Freq 2
n=0:1:600;                  % Sample vector
NFFT = 1024;                % FFT vector

% Discrete Time Domain Function
Xn = sin(2*pi*n*f1/fs) + sin(2*pi*n*f2/fs);

% Compute FFT and F-axis vector
fxn=fftshift(abs(fft(Xn,NFFT)));
if(rem(NFFT,2) == 0)
  % N is even
  f = linspace(-fs/2, (fs/2) - (fs/NFFT), NFFT);
else
  % N is odd
  f = linspace(-fs/2, (fs/2), NFFT);
end

%Plot Time Domain Representation
figure;
subplot(3,1,1);
plot(n,Xn);
grid on;
title(sprintf('x[n] = sin(2*pi*%d*n*Ts)+sin(2*pi*%d*n*Ts)',f1,f2),'fontsize',14);
xlabel('Time (Samples)','fontsize',14);
ylabel('Amplitude','fontsize',14);

subplot(3,1,2);
plot(f,fxn);
grid on;
title(sprintf('%d Point FFT: |X[n]|',NFFT),'fontsize',14);
xlabel('Time (Samples)','fontsize',14);
ylabel('Amplitude','fontsize',14);

subplot(3,1,3);
plot(f,20*log10(fxn));
grid on;
title(sprintf('%d Point FFT: |X[n]| - Log Scale',NFFT),'fontsize',14);
xlabel('Time (Samples)','fontsize',14);
ylabel('Amplitude (dB)','fontsize',14);

%%%%%%%%%%%%%%% PART 2 %%%%%%%%%%%%%%%
Y = conv(Xn,b);
figure;

% Compute FFT and F-axis vector
fyn=fftshift(abs(fft(Y,NFFT)));
if(rem(NFFT,2) == 0)
  % N is even
  f = linspace(-fs/2, (fs/2) - (fs/NFFT), NFFT);
else
  % N is odd
  f = linspace(-fs/2, (fs/2), NFFT);
end

% Show Visual
subplot(3,1,1);
plot(f,fxn,f,fftshift(abs(fft(250*b,NFFT))),'r');
grid on;
title('1024 point FFT: H[f] & X[f]');
xlabel('Frequency (Hz)');
ylabel('Amplitude');
legend('X[f]','H[f]')
subplot(3,1,1);

% Show Convolved overlay
subplot(3,1,2);
plot(f,fyn,f,fftshift(abs(fft(250*b,NFFT))),'r');
grid on;
title('1024 point FFT: h[n] * x[n] (with filter overlay)');
xlabel('Frequency (Hz)');
ylabel('Amplitude');
legend('X[f]','H[f]')
subplot(3,1,1);

% Show Convolved Ouput
subplot(3,1,3);
plot(f,fyn);
grid on;
% axis([-fs/2 fs/2 -80 60]);
title('1024 point FFT: h[n] * x[n]');
xlabel('Frequency (Hz)');
ylabel('Amplitude');
legend('X[f]','H[f]')
subplot(3,1,1);


